Numerical methodology to evaluate unipolar saturation current limit of DC corona discharge in complex geometries

The unipolar saturation current limit (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${I}_{sat}$$\end{document}Isat) gives an upper limit to the corona current that can be obtained from a unipolar corona discharge. Therefore, it implies a theoretical limit to the performance of unipolar corona discharge devices. However, it has not been widely used in practice because it is difficult to deal with complex discharge configurations in an analytical way. This study aims to establish and validate a numerical methodology to evaluate the maximum current, which numerically imitates the unipolar saturation current limit. It was shown that the maximum current has the same mathematical definition as the unipolar saturation current. For validation, the maximum current was compared with an analytical solution of the Poisson equation for the coaxial cylinders configuration. The differences between the maximum current and unipolar saturation current limit for the coaxial cylinders, pin-to-plane, and single wire-to-plane configurations were discussed in terms of the assumptions used in the semi-analytical derivation of the unipolar saturation current limit. The validated methodology was applied to a multiple wire-to-plane configuration, for which a semi-analytical expression of the unipolar saturation current limit has not yet been developed. The effects of geometric and operation parameters on the maximum currents of the multiple wire-to-plane configuration were analyzed. The results were regressed into a single formula.

that the discharge mode shifts to the bipolar regime, the uniformity and the homogeneity of the unipolar ion distribution are lost, degrading the performance of the devices. For instance, the collection efficiency of electrostatic precipitators is significantly reduced when the discharge mode shifts from the unipolar to the bipolar regime because the charges acquired by particulate matters are neutralized 10 . Kim and Hwang 11 showed that the discharge mode transition from the unipolar to the bipolar lowered the electric-to-kinetic energy conversion efficiency and the thrust performance of electrohydrodynamic flow generators. Recently, Wang et al. 12 utilized corona discharges to achieve a fast electrostatic printing of multifunctional materials to be attached on the human skin. A higher applied voltage may be preferred to achieve a faster printing, whereas it is essential to maintain the discharge in the unipolar regime for homogeneous and stable printing. Therefore, it is expected that the maximum performance of the unipolar corona discharge devices is closely related to the unipolar-tobipolar discharge mode transition, or the maximum amount of the unipolar ions that can be obtained from the discharge conditions of interest.
In this regard, Sigmond 8 semi-analytically derived an upper limit to the corona current that can be obtained from a unipolar corona discharge, which is called the unipolar saturation current limit ( I sat ). If the experimentally measured corona current ( I exp ) exceeds I sat , then the corona discharge of interest is, at least partially, bipolar. In this case, I exp does not result solely from the movement of the unipolar ions; at least a portion of I exp results from electrons, or secondary ionization in the drift region, or both of them 8 . In other words, I sat implies a theoretical limit to the performance of unipolar corona discharge devices.
To date, the evaluation of I sat has been based on the unipolar charge drift formula, which describes the timedependent change in unipolar charge density along the paths of ions 8 . Using this formula, the estimation of I sat reduces to the estimation of the lengths of electric field lines distributed in the inter-electrode space 13 . Many studies have analytically calculated the shape and length of electric field lines to obtain expressions for I sat in the discharge configurations of coaxial cylinders 8 , pin-to-plane 8,14 , and single wire-to-plane [14][15][16] . However, it is difficult to derive analytical solutions for electric field line distributions in complex configurations. As a result, I sat has not been widely used in the field of corona discharges in spite of its potential usefulness.
Recently, Kim et al. 17 suggested an approach to obtain I sat without using the unipolar charge drift formula, by proposing the concept of maximum current ( I max ). They described the calculation procedure to obtain I max and showed that I max had a similar value to I sat for a single wire-to-plane type corona discharge. Later, Kim and Hwang 11 simplified the calculation procedure by omitting the step of obtaining the size of the ionization region, which requires several iterative simulations.
This study aims to establish, validate, and apply a numerical methodology to evaluate I max for a complex discharge configuration. The calculation procedure to obtain I max was adopted from the method proposed by Kim and Hwang 11 . The procedure was validated by comparing the calculated values of I max with an analytical solution of the Poisson equation for coaxial cylinders configuration. After validation, this methodology was applied to pin-to-plane and single wire-to-plane configurations, and the values of I max were compared with those of I sat given by Sigmond 8,15 . Finally, I max of a multiple wire-to-plane configuration was evaluated. The dependence of I max on various geometric parameters and applied voltage was examined for this configuration.

Numerical methods
Unipolar corona discharge simulation. Corona discharge is a complex multi-physical phenomenon, and various processes are underlying, including ionization of neutral molecules by collision with electron, chemical reactions between numerous charged and neutral species, and transports of those discharge-induced species by electric field, ambient flow and ionic wind. Accordingly, numerical modeling of corona discharge in principle requires transport equations for electron, electron energy, and each species, and a set of reactions describing the chemical behavior of each species 5,18 . For a unipolar corona discharge, however, a drastically simplified model can be employed. The simplified model considers only one species: air ion. The air ion is a fictitious unipolar ion which represents various charged species generated from the unipolar corona discharge. Thus, all electrical behaviors of charged species in the unipolar corona discharge are implied in the behavior of the air ion. Within the simplified model, the generation and transport of neutral species are not considered. Therefore, chemical reactions occurring between various charged and neutral species are also not solved. Instead, the effect of the chemical reactions is implicitly included in the boundary condition of unipolar ion charge density. In the same manner, the behavior of electrons is also not considered. Thus, the potential effects of high gradient of electric field on the electron drift velocity and local field approximations cannot be discussed with the simplified model.
The simplified model consists of the steady-state governing equations for the electric field and distribution of the air ion given as follows: where V is the electric potential, q is the space charge density of the unipolar air ions generated by the corona discharge, − → j is the unipolar corona current density, µ is the electrical mobility of the ions, and − → E = −∇V is the electric field.
For the boundary condition of V , constant values were imposed to the surfaces of discharge ( V = V a ) and ground ( V = 0 ) electrodes, respectively. To specify the boundary condition of q at the discharge electrode, the ionization region was assumed to be circular and concentric with the center of curvature of the discharge electrode. The radius of the ionization region ( R i ) was assumed to be constant. For a pin-type discharge electrode, www.nature.com/scientificreports/ the ionization region was assumed to be a sphere with radius R i , concentric with the center of curvature of the tip of the pin. Similarly, for a wire-type discharge electrode, the ionization region was assumed to be a cylinder with radius R i , coaxial with the central axis of the wire. Then, a uniform and constant charge density q i was specified in the ionization region. At the ground electrode, the zero Neumann condition ( ∂q/∂n = 0 ) was imposed. After the V and q fields converged, the corona current was calculated by integrating − → j over the ground electrode.
Maximum current calculation. As mentioned in "Unipolar corona discharge simulation" section, in the simplified model, chemical aspects of corona discharge such as ionization are implicitly included in the boundary value of q at the discharge electrode, q i . Therefore, the problem of solving the unipolar corona discharge using the simplified model is equivalent to the problem of finding an appropriate value of q i ; the remaining is to solve the transport equations (Eqs. (1) and (2)) fluid-mechanically using the boundary condition. Accordingly, the unipolar corona current calculated by numerical simulation with the simplified model ( I num ) for a given applied voltage V a can be expressed as a function of q i as follows: If a value of q ′ i is appropriately determined to represent a real corona discharge, f is expected to be close to the experimentally measured corona current, that is, The unipolar saturation current I sat is an upper limit to the corona current obtainable from the unipolar corona discharge for a given V a . Intuitively, I num is expected to increase with increasing q i . Therefore, in the viewpoint of Eq. (3), I sat can be expressed as follows: Kim and Hwang 11 numerically solved the simplified model (Eqs. (1) and (2)) of the unipolar corona discharge. They fixed the value of R i and calculated I num by increasing q i stepwise by q i . Their results showed that when q i was sufficiently high, I num saturated at a certain value and did not increase even though q i increased further. Kim and Hwang 11 defined the saturated value of I num as the maximum current, I max . This description is exactly identical to Eq. (4). Therefore, I max can also be written as follows: In other words, I sat and I max are the semi-analytical and numerical implementations of the same expression, respectively. Thus, I sat and I max are expected to be equal to each other. In this study, I max was calculated by following the calculation procedure suggested by Kim and Hwang 11 . The flowchart of the calculation procedure is given in Supplementary Fig. S1. Sigmond 8,15 proposed the semi-analytical expressions of I sat for the coaxial cylinders (Eq. (6)), pin-to-plane (Eq. (7)), and single wire-to-plane (Eq. (8)) configurations as follows:

Analytical validation
Feng 19 reported a fully-analytical solution of the Poisson equation for corona discharge for the cylinders configuration using one-dimensional symmetricity. For other geometries, fully-analytical solutions are not available. Therefore, in this study, the calculated value of I max was validated using the current-voltage ( I-V ) relationship of the coaxial cylinders configuration. In "Analytical current-voltage curve for coaxial cylinders configuration" section, the fully-analytical I-V relationship of the coaxial cylinders configuration reported in Feng 19 was introduced. In "Maximum current calculation for coaxial cylinders configurations" section, I max was calculated for the coaxial cylinders configuration. After that, the calculated values of I max were compared with I sat and fullyanalytically computed corona current, which will be denoted by I analyt .  1) and (2)), the analytical I-V relationship for the coaxial cylinders configuration is given as follows 19 : www.nature.com/scientificreports/ where R w and R o are the radii of the inner wire and outer cylinder electrodes, respectively, and I analyt is the analytical corona current per unit length along the axial direction given as follows: C 1 is a non-negative integral constant of the Poisson equation and expressed as follows: where E w is the electric field strength at the surface of the inner wire electrode. Several studies employed the same approach with more details on the behaviors of electrons and ions such as ionization and attachment, and presented similar expressions to Eq. (9) 22,23 .
Because all other parameters are supposed to be given, Eq. (9) calculates the value of V a required to generate the corona current I analyt if the integral constant C 1 is given. To calculate C 1 , the value of E w should first be determined. The simplest and most common method to determine E w is to apply the Kaptzov hypothesis, which states that for all V a greater than the corona onset voltage, the value of E w remains constant at the onset value E on that is usually calculated by the Peek equation. The Peek equation for a cylindrical discharge electrode at the standard temperature and pressure is given as follows: where f is the roughness factor of the electrode surface.
Two coaxial cylinders configurations were adopted for validation: R w = 0.1 mm and R o = 25.5 mm 24 , and R w = 3.19 mm and R o = 290.5 mm 25 (see Fig. 1a). Figure 2 shows the comparisons between I exp and I analyt in the coaxial cylinders configurations. The value of µ was set to 2.0 × 10 -4 m 2 V −1 s −1 as reported 24,25 . The fully-analytical I-V curve given by Feng 19 matched with the experimental values 24,25 with f = 0.9 and 0.6, respectively. The value of f is 1.0 for the polished clean surface and 0.5-0.7 for dirty scratched wires 26 . Figure 2 also shows the I sat curves for each case. For the case of R w = 0.1 mm and R o = 25.5 mm, I exp was lower than I sat . It should be noted that I exp lower than I sat does not mean a unipolar corona discharge. For the case of R w = 3.19 mm and R o = 25.5 mm, I exp exceeded I sat at high V a . Therefore, it can be said that for the operation condition at which I exp was higher than I sat , at least a portion of I exp resulted from a bipolar corona discharge.
Maximum current calculation for coaxial cylinders configurations. Following the procedure described in "Maximum current calculation" section, I max was calculated for the coaxial cylinders configurations. Figure 3a shows the saturation profiles of I num versus q i for positive discharge under the conditions R w = 0.1 mm, R o = 25.5 mm, and R i = 3 R w . For any value of V a , I num became saturated as q i increased (see square markers in Fig. 3a), as in Kim and Hwang 11 . Thus, the value of I num corresponding to the square marker can be accepted as I max for a given V a . Figure 3a also shows (right vertical axis) the profiles of the reduced electric field strength at the surface of the discharge electrode ( E w /N , where N is the number density of neutral gas molecules) versus q i . When q i increased, E w /N decreased while I num increased. At any saturation point, E w /N was below 120 Td (1 Td = 10 −21 V m 2 ), which is equivalent to the breakdown electric field strength ( E b ) of air at the standard temperature and pressure. The ionization coefficient is greater than the attachment coefficient where 27,28 . Electrons are sufficiently accelerated and energetic enough to ionize the colliding neutral molecules when E > E b . On the other hand, when E < E b , electrons are not sufficiently accelerated. Insufficiently energetic electrons cannot ionize neutral molecules, but they are attached to the molecules. Therefore, it is reasonable to assume that E w should also remain greater than E b even if the unipolar corona discharge is in its saturated state. The circular markers in Fig. 3a indicate E w /N = 120 Td, and the cross markers indicate the corresponding I num . Hereinafter, if E w /N reached 120 Td, the corresponding I num was accepted as I max even though I num was not yet completely saturated.
As mentioned in "Analytical current-voltage curve for coaxial cylinders configuration" section, the value of E w should be first specified to calculate I analyt . As discussed with Fig. 2, E on was used as E w for a real corona discharge. However, for the saturated state, the representative value of E w has not been reported. In other words, it is unclear which value should be used as E w in Eq. (11) to calculate I analyt representing I sat and I max . Obviously, the Kaptzov hypothesis cannot be employed for the saturated state, as shown in Fig. 2; I analyt calculated by using E w = E on was significantly different from I sat for the case R w = 0.1 mm and R o = 25.5 mm. Sigmond 8 derived the semi-analytical expression of I sat without using information on E w . Instead, an average electric field strength (e.g., E = V a /R o for the coaxial cylinders configuration) was used. Therefore, without a better option, the breakdown electric field strength E b was used as E w for the saturated state, as also used for I max . I analyt,b given in Fig. 3b indicate the value of I analyt calculated using E w = E b . As shown in Fig. 3b, I max was the same as I analyt,b for all V a .
The dashed lines in Fig. 3b indicate the I sat curves. I max and I sat had the similar magnitudes as in Kim et al. 17 , but with some deviations. The semi-analytical expressions for I sat given in Eqs. (6)-(8) were originally derived under various assumptions such as the Deutsch assumption, which states that the space charge does not distort www.nature.com/scientificreports/ the shape of the space-charge-free (i.e. Laplacian) electric field lines, but only scale the strength of the electric field; and the small discharge electrode assumption, which can be expressed as R w /R o ≪ 1 8 . Because the Deutsch assumption is valid for the coaxial cylinders configuration, the deviation between I sat and I max might result from the invalidity of the small discharge electrode assumption. For simple comparison, I max was assumed to have the same form as I sat , i.e. I max = C 2 V γ a , where C 2 is an arbitrary constant. Then, the exponent γ was obtained  In other words, the smaller R w /R o was, the closer γ was to the semi-analytical value, 2 (Eq. (6)), which was consistent with the small discharge electrode assumption. The effects of three different R i values (1.5 R w , 3 R w , and 5 R w ) on I max were tested. For all conditions, R i had a negligible effect on I max , in agreement with Kim and Hwang 11 . For a real corona discharge, the shape and size of the ionization region are probably not constant, and the distribution of charges in the ionization region is also not uniform. Thus, numerical modeling of the real corona discharge should consider those effects. However, I max is an attempt to find a limit value of corona current obtainable under a certain applied voltage from a certain discharge configuration, not a real physical current. In the viewpoint of Eq. (5), I max is a current calculated with an infinitely high q i . When the value of q i is infinitely high, the way in which the boundary condition of q i is imposed should not affect the final value of I max .

Results and discussion
The validated calculation procedure was applied to the pin-to-plane and single wire-to-plane configurations to calculate I max . The results were compared with I sat curves. Then, the same methodology was applied to the multiple wire-to-plane configuration to examine the effects of various geometric parameters and applied voltage on I max . The value of µ was set to 1.4 × 10 -4 and 1.6 × 10 -4 m 2 V -1 s -1 for positive and negative discharges, respectively 29 .
Pin-to-plane configuration. The computational domain for the positive corona discharge in the pinto-plane configuration was adopted from Adamiak and Atten 30 (see Fig. 1b). Figure 4a shows that the values of I max were approximately the same as those of I sat in various pin-to-plane configurations. With assuming  www.nature.com/scientificreports/ I max = C 2 V γ a as in Section Maximum current calculation for coaxial cylinders configurations, the values of γ were given by 2.047, 2.089, and 2.090 for R t /L = 0.0011, 0.0031, and 0.0048, respectively, which were very close to the semi-analytical value, 2 (Eq. (7)). As mentioned in Section Analyticalcurrent voltage curve for coaxialcylinders configuration Based and Maximum current calculation for coaxial cylinders configurations, the semianalytical expressions of I sat given in Eqs. (6)- (8) were derived under the small discharge electrode assumption. Therefore, it can be said that the pin electrodes with R t /L = 0.0011 to 0.0048 used in this study were small enough to neglect the effect of R t on I max . It is noteworthy that unlike the coaxial cylinders configurations, the Deutsch assumption is basically not valid for the non-symmetric geometry such as the pin-to-plane and wire-to-plane configurations 31,32 . Jones and Davies 31 concluded that the use of the Deutsch assumption to the pin-to-plane configuration should be avoided because it can cause a significant error. Nevertheless, I max was in good agreement with I sat for the considered conditions. It is unclear and has not been reported in the literature whether and how the validity of the Deutsch assumption affects I sat . However, the analysis of the theoretical background of I sat is not in the scope of this study, and thus, will not be further discussed. The change in R i (1.5 R t , 3 R t , and 5 R t ) had little effect on I max of the pin-to-plane configuration as shown in Supplementary Fig. S2.
Single wire-to-plane configuration. The computational domain for the negative corona discharge in the single wire-to-plane configuration was adopted from Ohkubo et al. 33 (see Fig. 1c). Figure 4b shows I max versus V a under D = 100 mm and various R w . The values of γ were greater than 2 for all R w and increased to greater than 3 as R w increased. Because corona discharge would be less intensive for thicker wires, I max is also expected to be smaller for thicker wires. In addition, this result was consistent with the small discharge electrode assumption; the smaller R w was, the closer I max was to I sat . As in other configurations, R i had little effect on I max of the single wire-to-plane configuration (see Supplementary Fig. S3). It is noted that both I sat given in Eq. (8) and I max discussed within the single wire-to-plane configuration are the current values for one plane.
Multiple wire-to-plane configuration. Figure 1d shows the multiple wire-to-plane configuration. D and S are the wire-to-ground distance and half of the wire-to-wire spacing, respectively. The normalized position along the ground electrode is represented by x/S , and θ is the angle of incidence. To observe the effects of various combinations of S and D , the following parameter was used: As mentioned in Section Single wire to plane configuration, I sat = 1.62µεV 2 a /D 2 for the single wire-to-plane configuration given in Eq. (8) describes the saturation current per one plane. For the multiple wire-to-plane configuration, this is equivalent to one cell (green shading in Fig. 1d). Assuming that I max in the multiple wireto-plane configuration exhibits similar behavior to that of the single wire-to-plane configuration, the expression of I max for one cell of the multiple wire-to-plane configuration will have the similar form with that for the single wire-to-plane, that is, I max,cell = C 3 µεV γ a /D β , where C 3 , γ , and β are constants. Then, I max of the entire multiple wire-to-plane configuration can be expressed as I max,multiple = 2n w I max,cell = 2n w C 3 µεV γ a /D β , where n w is the number of wires. In addition, the wire-to-wire spacing S would affect the magnitude of I max by an electrical shielding. Therefore, I max,multiple will have the following form: where g(S) is a correlation function that describes the dependence of I max on S . The actual simulation was conducted over half of one cell, or a unit (blue shading in Fig. 1d), using a symmetry boundary condition. I max for one unit of the multiple wire-to-plane configuration ( I max,unit ) was obtained by integrating the current density over S (see Fig. 1d). Then, I max for one cell was calculated by I max,cell = 2I max,unit . Hereinafter, I max for the multiple wire-to-plane configuration refers to I max,cell to directly compare the maximum current with the saturation current. Thus, I max for the multiple wire-to-plane configuration is expressed as follows: where C 4 is a constant. The ranges of S , D , and V a considered were 50-1000 mm, 50-250 mm, and 45-70 kV, respectively. For simplicity, the dependence of I max on R w was not considered, and the value of R w was fixed at 1.59 mm. Figure 5a shows I max calculated for various values of S and V a with D = 114.3 mm. When S was small so that θ S/D was less than approximately 65°, I max increased with increasing S . However, as S increased further so that θ S/D became greater than 65°, I max converged to I sat of the single wire-to-plane configuration. When θ S/D > 65°, two adjacent wires were too far apart to affect each other. They then behaved as two independent single wire-to-plane configurations. For D = 114.3 mm, S = 250, 500, and 1000 mm gave θ S/D = 65.4, 77.1, and 83.5°, respectively. In other words, for D = 114.3 mm, the multiple wire-to-plane configurations with S = 250, 500, and 1000 mm were equivalent to simple series of n w single wire-to-plane configurations. Figure 5a also shows I sat of the single wire-to-plane configuration for D = 114.3 mm. For V a = 70 kV, I max for the multiple wire-to-plane configuration converged closely to I sat for the single wire-to-plane configuration as S increased. On the other hand, for lower V a , the converged value of I max was smaller than I sat . This is consistent with the results for the single wire-to-plane configuration with thick wires demonstrated in Fig. 4b.
The magnitude of the corona current density ( j ) on the ground electrode is known to follow the Warburg distribution, which is given as follows 8 : www.nature.com/scientificreports/ where θ = tan −1 (x/D) is the angle of incidence (see Fig. 1d), and m is an empirical constant. Kim et al. 17 stated that for the single wire-to-plane configuration, the distributions of j corresponding to I exp and I num , namely j exp and j num , followed the Warburg distribution with m = 4.2. Figure 5b shows that for S larger than 250 mm, the distribution of j corresponding to I max (that is, j max ) also followed the Warburg distribution with m = 4.2. Figure 5b also shows that for S larger than 250 mm, j max dropped to zero before reaching the midpoint between two wires ( x/S = 1). In this situation, the two wires did not influence each other, and each cell in the multiple wire-to-plane configuration converged to the single wire-to-plane configuration, as shown in Fig. 5a. For the single wire-to-plane configuration, ions generated from a discharge electrode do not drift outside θ > 65°; hence, j drops to zero at approximately θ = 65°8 ,32 . For S less than 250 mm, however, j max did not fall to zero before reaching x/S = 1 because of the electrical shielding by the neighboring wires 34 . Therefore, in this range, I max was smaller by the undrawn area of the j max distribution, resulting in a decrease in I max with decreasing S , as shown in Fig. 5a. In addition, for S less than 250 mm, j max itself deviated from the Warburg distribution. When S = 50 mm, j max at x/S = 1 was approximately 33% smaller than its value in the Warburg distribution. This further reduced the value of I max . Figure 6 shows the values of I max normalized by the value of I max at S = 1000 mm for various V a . Note that the values of I max calculated for S = 1000 mm approximate the values of I max of the single wire-to-plane configuration for the considered range of D (50-250 mm). Figure 5a showed that the magnitude of I max increased with increasing V a . However, Fig. 6 revealed that the normalized magnitude of I max did not depend on V a ; I max for different V a collapsed into a single profile. For θ S/D < 45°, the normalized I max was linearly proportional to θ S/D . Regression analysis in the linear region (Regime I) yielded I max /I max (S = 1000) = 1.18θ S/D − 0.065 , where θ S/D is in radians (solid line in Fig. 6). For θ S/D > 65° (Regime III), I max was the same as that in the single wire-to-plane configuration, as also shown in Figs. 5a and b. Then, replacing I max (S = 1000) using Eq. (15), the expression for I max is given by where θ S/D is in radians in the regression equation. Regime II (45° < θ S/D < 65°) represents a transition between Regimes I and III. Using regression, the three regimes can also be represented in a continuous manner as follows (dashed curve in Fig. 6): where θ S/D is in radians. Figure 7 shows the dependence of I max on D for various values of S and V a . Taking the logarithm with respect to D of both sides of Eq. (15) gave the slope of the least-squares fit line as β . For a given S , β was approximately independent of V a . The average values of β were 3.06, 2.76, 2.66, 2.54, 2.40, 2.34, and 2.24 for S = 50, 80, 114.3, 150, 200, 250, and 500 mm, respectively. That is, the dependence of I max on D was affected by the value of S ; the value of β decreased and converged to that of the single wire-to-plane configuration as S increased.
Taking the logarithm with respect to V a of both sides of Eq. (15) gave the slope of the least-squares fit line as γ . Figure 8 shows a plot of γ versus θ S/D for given D . Note that for each value of D , an increase in θ S/D corresponds www.nature.com/scientificreports/ to an increase in S . Overall, the value of γ ranged in approximately 2.5-3.0 and decreased slightly with increasing θ S/D or S.

Conclusion
In this study, a numerical methodology to evaluate the unipolar saturation current of DC corona discharges was proposed using the concept of maximum current. It was shown that the maximum current has the same mathematical definition as the unipolar saturation current. For validation, the maximum current of the unipolar corona discharge for the coaxial cylinders configuration was compared with the analytical solution of the Poisson equation. For the coaxial cylinders and single wire-to-plane configurations, the difference between the maximum current and the unipolar saturation current decreased as the radius of the wire electrode decreased,  www.nature.com/scientificreports/ confirming the effect of small discharge electrode assumption. For the pin-to-plane configuration, the maximum current was approximately the same as the unipolar saturation current. The calculation procedure was applied to a multiple wire-to-plane configuration, for which a semi-analytical expression of the saturation current has not yet been reported. The results showed that the maximum current was influenced by the electrical shielding from the neighboring discharge electrodes. The effects of geometric and operation parameters on the maximum currents of the multiple wire-to-plane configuration were examined, showing a potential usage for complex discharge configurations.

Data availability
The datasets used and/or analyzed during the current study available from the corresponding author on reasonable request.